home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-04 / amsf20.zip / DGELS.FOR < prev    next >
Text File  |  1992-01-06  |  6KB  |  187 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DGELS
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH
  8. C           SYMMETRIC COEFFICIENT MATRIX UPPER TRIANGULAR PART OF WHICH
  9. C           IS ASSUMED TO BE STORED COLUMNWISE.
  10. C
  11. C        USAGE
  12. C           CALL DGELS(R,A,M,N,EPS,IER,AUX)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           R      - DOUBLE PRECISION M BY N RIGHT HAND SIDE MATRIX
  16. C                    (DESTROYED). ON RETURN R CONTAINS THE SOLUTION OF
  17. C                    THE EQUATIONS.
  18. C           A      - UPPER TRIANGULAR PART OF THE SYMMETRIC DOUBLE
  19. C                    PRECISION M BY M COEFFICIENT MATRIX.  (DESTROYED)
  20. C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
  21. C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
  22. C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS
  23. C                    RELATIVE TOLERANCE FOR TEST ON LOSS OF
  24. C                    SIGNIFICANCE.
  25. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  26. C                    IER=0  - NO ERROR,
  27. C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
  28. C                             PIVOT ELEMENT AT ANY ELIMINATION STEP
  29. C                             EQUAL TO 0,
  30. C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
  31. C                             CANCE INDICATED AT ELIMINATION STEP K+1,
  32. C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
  33. C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
  34. C                             ABSOLUTELY GREATEST MAIN DIAGONAL
  35. C                             ELEMENT OF MATRIX A.
  36. C           AUX    - DOUBLE PRECISION AUXILIARY STORAGE ARRAY
  37. C                    WITH DIMENSION M-1.
  38. C
  39. C        REMARKS
  40. C           UPPER TRIANGULAR PART OF MATRIX A IS ASSUMED TO BE STORED
  41. C           COLUMNWISE IN M*(M+1)/2 SUCCESSIVE STORAGE LOCATIONS, RIGHT
  42. C           HAND SIDE MATRIX R COLUMNWISE IN N*M SUCCESSIVE STORAGE
  43. C           LOCATIONS. ON RETURN SOLUTION MATRIX R IS STORED COLUMNWISE
  44. C           TOO.
  45. C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
  46. C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
  47. C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
  48. C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
  49. C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
  50. C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
  51. C           GIVEN IN CASE M=1.
  52. C           ERROR PARAMETER IER=-1 DOES NOT NECESSARILY MEAN THAT
  53. C           MATRIX A IS SINGULAR, AS ONLY MAIN DIAGONAL ELEMENTS
  54. C           ARE USED AS PIVOT ELEMENTS. POSSIBLY SUBROUTINE DGELG (WHICH
  55. C           WORKS WITH TOTAL PIVOTING) WOULD BE ABLE TO FIND A SOLUTION.
  56. C
  57. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  58. C           NONE
  59. C
  60. C        METHOD
  61. C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
  62. C           PIVOTING IN MAIN DIAGONAL, IN ORDER TO PRESERVE
  63. C           SYMMETRY IN REMAINING COEFFICIENT MATRICES.
  64. C
  65. C     ..................................................................
  66. C
  67.       SUBROUTINE DGELS(R,A,M,N,EPS,IER,AUX)
  68. C
  69. C
  70.       DIMENSION A(1),R(1),AUX(1)
  71.       DOUBLE PRECISION R,A,AUX,PIV,TB,TOL,PIVI
  72.       IF(M)24,24,1
  73. C
  74. C     SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT
  75.     1 IER=0
  76.       PIV=0.D0
  77.       L=0
  78.       DO 3 K=1,M
  79.       L=L+K
  80.       TB=DABS(A(L))
  81.       IF(TB-PIV)3,3,2
  82.     2 PIV=TB
  83.       I=L
  84.       J=K
  85.     3 CONTINUE
  86.       TOL=EPS*PIV
  87. C     MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.
  88. C     PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
  89. C
  90. C
  91. C     START ELIMINATION LOOP
  92.       LST=0
  93.       NM=N*M
  94.       LEND=M-1
  95.       DO 18 K=1,M
  96. C
  97. C     TEST ON USEFULNESS OF SYMMETRIC ALGORITHM
  98.       IF(PIV)24,24,4
  99.     4 IF(IER)7,5,7
  100.     5 IF(PIV-TOL)6,6,7
  101.     6 IER=K-1
  102.     7 LT=J-K
  103.       LST=LST+K
  104. C
  105. C     PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
  106.       PIVI=1.D0/A(I)
  107.       DO 8 L=K,NM,M
  108.       LL=L+LT
  109.       TB=PIVI*R(LL)
  110.       R(LL)=R(L)
  111.     8 R(L)=TB
  112. C
  113. C     IS ELIMINATION TERMINATED
  114.       IF(K-M)9,19,19
  115. C
  116. C     ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
  117. C     ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.
  118.     9 LR=LST+(LT*(K+J-1))/2
  119.       LL=LR
  120.       L=LST
  121.       DO 14 II=K,LEND
  122.       L=L+II
  123.       LL=LL+1
  124.       IF(L-LR)12,10,11
  125.    10 A(LL)=A(LST)
  126.       TB=A(L)
  127.       GO TO 13
  128.    11 LL=L+LT
  129.    12 TB=A(LL)
  130.       A(LL)=A(L)
  131.    13 AUX(II)=TB
  132.    14 A(L)=PIVI*TB
  133. C
  134. C     SAVE COLUMN INTERCHANGE INFORMATION
  135.       A(LST)=LT
  136. C
  137. C     ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT
  138.       PIV=0.D0
  139.       LLST=LST
  140.       LT=0
  141.       DO 18 II=K,LEND
  142.       PIVI=-AUX(II)
  143.       LL=LLST
  144.       LT=LT+1
  145.       DO 15 LLD=II,LEND
  146.       LL=LL+LLD
  147.       L=LL+LT
  148.    15 A(L)=A(L)+PIVI*A(LL)
  149.       LLST=LLST+II
  150.       LR=LLST+LT
  151.       TB=DABS(A(LR))
  152.       IF(TB-PIV)17,17,16
  153.    16 PIV=TB
  154.       I=LR
  155.       J=II+1
  156.    17 DO 18 LR=K,NM,M
  157.       LL=LR+LT
  158.    18 R(LL)=R(LL)+PIVI*R(LR)
  159. C     END OF ELIMINATION LOOP
  160. C
  161. C
  162. C     BACK SUBSTITUTION AND BACK INTERCHANGE
  163.    19 IF(LEND)24,23,20
  164.    20 II=M
  165.       DO 22 I=2,M
  166.       LST=LST-II
  167.       II=II-1
  168.       L=A(LST)+.5D0
  169.       DO 22 J=II,NM,M
  170.       TB=R(J)
  171.       LL=J
  172.       K=LST
  173.       DO 21 LT=II,LEND
  174.       LL=LL+1
  175.       K=K+LT
  176.    21 TB=TB-A(K)*R(LL)
  177.       K=J+L
  178.       R(J)=R(K)
  179.    22 R(K)=TB
  180.    23 RETURN
  181. C
  182. C
  183. C     ERROR RETURN
  184.    24 IER=-1
  185.       RETURN
  186.       END
  187.